A Revised Abaqus® Procedure for Fracture Path Simulation Based on the Material Effort Criterion

This paper presents the results of computer simulations of fracture in three laboratory tests: the three-point bending of a notched beam cut from sandstone, the pull-out test of a self-undercutting anchor fixed in sandstone, and the pull-out test of a bar embedded in concrete. Five material failure criteria were used: Rankine, Coulomb–Mohr, Drucker–Prager, Ottosen–Podgórski, and Hoek–Brown. These criteria were implemented in the Abaqus® FEA system to work with the crack propagation modeling method—extended finite element method (X-FEM). All criteria yielded similar force–displacement relationships and similar crack path shapes. The improved procedure gives significantly better, close-to-real crack propagation paths than can be obtained using the standard subroutines built into the Abaqus® system.


Introduction
Numerical fracture simulation is a convenient and frequently used method for determining the critical states of materials and structures.However, it requires specialized software that, in addition to determining the critical load values, can determine the shape and range of the propagating crack.This is especially important when studying brittle or quasi-brittle materials, where the singularities of the stress field near the crack tip are the cause of large errors in determining the parameters of the propagation process.Such materials can include, for example, concrete and hard rocks, such as sandstone, granite, and dolomite, for which the fracture process is characterized by great dynamics.These are materials often tested in pull-out tests to determine the safety and load-bearing capacity of anchors fixed in walls or the size of a chipped rock fragment and the force required to break it off [1][2][3][4][5].This is one of the possible rock removal methods tested for mine rescue [6].One of the software tools that can be used for these simulations is the advanced and popular Simulia Abaqus ® system [7], which offers techniques to perform this type of study.Using this system in a 'standard' way, however, comes with some limitations that, under certain conditions, produce simulation results that are not very accurate in relation to predicting the direction of crack propagation.These limitations can be overcome thanks to the openness of the Abaqus ® system, which enables the implementation of independently developed subroutines that control the calculations.
This paper continues the authors' work [8] on their own method for simulating crack propagation in the Abaqus ® FEA system.This method uses selected material failure criteria and determines the crack initiation and propagation angle.It cooperates with the X-FEM method to simulate the crack propagation.The authors' method was implemented in Fortran as an Abaqus ® User Subroutine [7].The main goal of this research is to compare several material failure criteria that were implemented in the Abaqus ® system by the authors and to check if it is possible to easily implement any criterion to work with the Abaqus ® system.The X-FEM method is one of the most popular fracture simulation methods used in FEA systems [9].The crack is implemented by finite element separation using an enrichment function added to the element shape function, which allows the simulation of displacement discontinuities in the element area.Typically, programs support a single crack jumping from edge to edge of an element during simulation in two-dimensional models, but solutions to three-dimensional problems can be found in the literature [10] and even provide crack branching [11].
Simulia's Abaqus ® FEA system enables complex simulations of nonlinear processes occurring in brittle and quasi-brittle materials during cracking.This program allows the use of many different crack simulation methods, including the X-FEM method in combination with the Cohesive Zone Method (CZM), which makes it possible to achieve a high degree of independence of the simulation results obtained from the FEM mesh density.The (CZM) technique of avoiding the singularity of the stress field in the vicinity of the crack tip, proposed by Barenblatt [12], is still frequently used in fracture simulation analyses of brittle materials and continues to be improved and developed by today's researchers [13,14].
The Abaqus ® system allows the selection of several default criteria for determining the crack initiation and propagation angle.The most typical criterion for applications in homogeneous materials, in which the crack can run in any direction, is the MAXPS criterion, i.e., the criterion of maximum principal stresses.There are several different criteria that cooperate with X-FEM formulation in Abaqus ® , such as MAXS (maximum nominal stress) or QUADS (quadratic nominal stress), but they are not fit for analyses of brittle materials-by selecting these criteria, the slot can only be routed in directions along the axes of the local finite elements.For the purposes of this paper, the MAXPS criterion will be called the default or built-in criterion.It works in the following way: the program reads the stresses at the integration points of the element to which the crack tip belongs, and the fracture propagation angle is equal to the average of the rotation angles of the stress tensors to the principal stresses at four integration points.Fracture occurs when the average maximum principal stress exceeds the specified tensile strength.A diagram of how this method works is shown in Figure 1.
Materials 2024, 17, x FOR PEER REVIEW 2 of 19 The X-FEM method is one of the most popular fracture simulation methods used in FEA systems [9].The crack is implemented by finite element separation using an enrichment function added to the element shape function, which allows the simulation of displacement discontinuities in the element area.Typically, programs support a single crack jumping from edge to edge of an element during simulation in two-dimensional models, but solutions to three-dimensional problems can be found in the literature [10] and even provide crack branching [11].
Simulia's Abaqus ® FEA system enables complex simulations of nonlinear processes occurring in brittle and quasi-brittle materials during cracking.This program allows the use of many different crack simulation methods, including the X-FEM method in combination with the Cohesive Zone Method (CZM), which makes it possible to achieve a high degree of independence of the simulation results obtained from the FEM mesh density.The (CZM) technique of avoiding the singularity of the stress field in the vicinity of the crack tip, proposed by Barenblatt [12], is still frequently used in fracture simulation analyses of brittle materials and continues to be improved and developed by today's researchers [13,14].
The Abaqus ® system allows the selection of several default criteria for determining the crack initiation and propagation angle.The most typical criterion for applications in homogeneous materials, in which the crack can run in any direction, is the MAXPS criterion, i.e., the criterion of maximum principal stresses.There are several different criteria that cooperate with X-FEM formulation in Abaqus ® , such as MAXS (maximum nominal stress) or QUADS (quadratic nominal stress), but they are not fit for analyses of brittle materials-by selecting these criteria, the slot can only be routed in directions along the axes of the local finite elements.For the purposes of this paper, the MAXPS criterion will be called the default or built-in criterion.It works in the following way: the program reads the stresses at the integration points of the element to which the crack tip belongs, and the fracture propagation angle is equal to the average of the rotation angles of the stress tensors to the principal stresses at four integration points.Fracture occurs when the average maximum principal stress exceeds the specified tensile strength.A diagram of how this method works is shown in Figure 1.The X-FEM technique implemented in Abaqus ® is quite immune to the density of the finite element mesh.Abaqus ® does not lead the stresses to infinity at the crack tip-the maximum value is close to the tensile strength; therefore, regardless of the size of the finite elements, the destructive forces (i.e., peak forces) in the model are very close to each other, which was proven in several previous author's works.However, the difference in mesh size may lead to differences in the crack line, as Abaqus ® always leads the crack from edge to edge of the element.This may cause a crack in a model with a very small number of finite elements that are not as smooth as in a model with a dense mesh.However, the The X-FEM technique implemented in Abaqus ® is quite immune to the density of the finite element mesh.Abaqus ® does not lead the stresses to infinity at the crack tip-the maximum value is close to the tensile strength; therefore, regardless of the size of the finite elements, the destructive forces (i.e., peak forces) in the model are very close to each other, which was proven in several previous author's works.However, the difference in mesh size may lead to differences in the crack line, as Abaqus ® always leads the crack from edge to edge of the element.This may cause a crack in a model with a very small number of finite elements that are not as smooth as in a model with a dense mesh.However, the simulations analyzed below have enough finite elements to consider the resulting cracks smooth enough.

Own Method for Selecting the Direction of Fracture Propagation
The method for predicting the direction of crack propagation was based on the authors' proposal using the direction of the material effort field gradient [15].It provides results almost identical (Figure 2) to the classical maximum energy release rate method [16] and is much more convenient and easier to apply for complex analyses [8].simulations analyzed below have enough finite elements to consider the resulting cracks smooth enough.

Own Method for Selecting the Direction of Fracture Propagation
The method for predicting the direction of crack propagation was based on the authors' proposal using the direction of the material effort field gradient [15].It provides results almost identical (Figure 2) to the classical maximum energy release rate method [16] and is much more convenient and easier to apply for complex analyses [8].This method for selecting the direction of crack propagation was implemented in the Abaqus ® system using the Abaqus ® User Subroutine and written in Fortran programming language.The most important subroutine used in the simulations is the UDMGINI (User Damage Initiation), which Abaqus ® runs on every iteration of the calculations (even several thousand times during the simulation).The input data of the subroutine are the stresses and coordinates at the integration points of the enriched elements and other available data that can be imported from other subroutines.The output data is only the value of the material effort and the crack propagation angle.If the subroutine determines the material effort to be greater than 1, Abaqus ® leads the crack to the next finite element.The authors also wrote two other auxiliary subroutines: UVARM (user subroutine to generate element output) and URDFIL (user subroutine to read the result file).
The subroutine codes and the operation of the authors' method for determining the direction of crack propagation are described in the authors' previous publication [8].The idea of how the algorithm works is simple.For each load increment, the program selects several dozen integration points around the crack tip, reads the stresses at these points, and then uses the stresses to calculate the material effort, which depends on the selected material failure criterion.The effort value is then interpolated to points lying equidistant from the crack tip. Figure 3b shows an example of a graph of the relationship between the values of the material effort at the integration points and the angle θ between the horizontal axis and the line connecting the crack tip with each integration point.The This method for selecting the direction of crack propagation was implemented in the Abaqus ® system using the Abaqus ® User Subroutine and written in Fortran programming language.The most important subroutine used in the simulations is the UDMGINI (User Damage Initiation), which Abaqus ® runs on every iteration of the calculations (even several thousand times during the simulation).The input data of the subroutine are the stresses and coordinates at the integration points of the enriched elements and other available data that can be imported from other subroutines.The output data is only the value of the material effort and the crack propagation angle.If the subroutine determines the material effort to be greater than 1, Abaqus ® leads the crack to the next finite element.The authors also wrote two other auxiliary subroutines: UVARM (user subroutine to generate element output) and URDFIL (user subroutine to read the result file).
The subroutine codes and the operation of the authors' method for determining the direction of crack propagation are described in the authors' previous publication [8].The idea of how the algorithm works is simple.For each load increment, the program selects several dozen integration points around the crack tip, reads the stresses at these points, and then uses the stresses to calculate the material effort, which depends on the selected material failure criterion.The effort value is then interpolated to points lying equidistant from the crack tip. Figure 3b shows an example of a graph of the relationship between the values of the material effort at the integration points and the angle θ between the horizontal axis and the line connecting the crack tip with each integration point.The interpolation polynomial is then fitted to points with these coordinates, and then the subroutine finds its local minimum.The angle for which the local minimum is obtained is the angle of the next step of crack propagation.Additionally, the crack propagates when the material effort at the integration point closest to the crack tip is greater than 1.This is a better method than the built-in method, where the crack is determined by the average from the four integration points.interpolation polynomial is then fitted to points with these coordinates, and then the subroutine finds its local minimum.The angle for which the local minimum is obtained is the angle of the next step of crack propagation.Additionally, the crack propagates when the material effort at the integration point closest to the crack tip is greater than 1.This is a better method than the built-in method, where the crack is determined by the average from the four integration points.
(a) (b) The method of selecting integration points taken into account for the polynomial fitting stage was developed by trial and error.There should be enough points to fit the polynomial as accurately as possible.These points should not be too close to the crack tip, as there are quite large disturbances.However, they should be close enough due to the fracture being determined by the stresses in a small area around the crack tip.It was decided that the algorithm would select the radius r so that there are 50 to 100 integration points at a distance from r to 1.5 r from the crack tip.Too large a number would also lead to an unnecessary increase in the duration of the simulation.Moreover, only integration points that lay a maximum of 120° to the left and right of the angle of the existing crack are taken into account in the calculations.This angle could not be too large to reject the stresses at points close to the already existing crack.It also could not be too small, or it would not provide enough integration points for calculations.An illustration of the method of selecting integration points is shown in Figure 3a.Also, by trial and error, it was decided that the polynomial fitting to the points would be of the fourth degree.Too low a degree of the polynomial caused the function to fit the points inaccurately in some cases.
The selection of integration points was made as best as possible, but some problems still arose.Figure 4a shows a crack running perpendicular to the edge of the model.The area of collecting integration points overlapped the model, so the fit of the 4th-degree polynomial to the material effort values was not problematic, which is also visible in Figure 4b.Unfortunately, when the crack approached the edge, the area of collecting integration points partially went outside the model (Figure 5a).For this reason, for a large range of the angle θ, there were no data of stresses at the integration points, as can be seen in Figure 5b.In this case, fitting the 4th-degree polynomial was incorrect, as can be seen in the figure.The local minimum suddenly jumped to the wrong angle, causing the crack in the simulations to turn incorrectly to the left or right.This problem was solved in such a way that when the algorithm encounters a situation where a large part of the area taken into account for collecting integration points goes beyond the model, the polynomial is changed to the 2nd degree.Figure 5b shows that the local minimum for the 2nd-degree polynomial was correct for an angle of −90°, and the crack still went vertically.The method of selecting integration points taken into account for the polynomial fitting stage was developed by trial and error.There should be enough points to fit the polynomial as accurately as possible.These points should not be too close to the crack tip, as there are quite large disturbances.However, they should be close enough due to the fracture being determined by the stresses in a small area around the crack tip.It was decided that the algorithm would select the radius r so that there are 50 to 100 integration points at a distance from r to 1.5 r from the crack tip.Too large a number would also lead to an unnecessary increase in the duration of the simulation.Moreover, only integration points that lay a maximum of 120 • to the left and right of the angle of the existing crack are taken into account in the calculations.This angle could not be too large to reject the stresses at points close to the already existing crack.It also could not be too small, or it would not provide enough integration points for calculations.An illustration of the method of selecting integration points is shown in Figure 3a.Also, by trial and error, it was decided that the polynomial fitting to the points would be of the fourth degree.Too low a degree of the polynomial caused the function to fit the points inaccurately in some cases.
The selection of integration points was made as best as possible, but some problems still arose.Figure 4a shows a crack running perpendicular to the edge of the model.The area of collecting integration points overlapped the model, so the fit of the 4th-degree polynomial to the material effort values was not problematic, which is also visible in Figure 4b.Unfortunately, when the crack approached the edge, the area of collecting integration points partially went outside the model (Figure 5a).For this reason, for a large range of the angle θ, there were no data of stresses at the integration points, as can be seen in Figure 5b.In this case, fitting the 4th-degree polynomial was incorrect, as can be seen in the figure.The local minimum suddenly jumped to the wrong angle, causing the crack in the simulations to turn incorrectly to the left or right.This problem was solved in such a way that when the algorithm encounters a situation where a large part of the area taken into account for collecting integration points goes beyond the model, the polynomial is changed to the 2nd degree.Figure 5b shows that the local minimum for the 2nd-degree polynomial was correct for an angle of −90 • , and the crack still went vertically.
There is also another problem that seems to have no simple solution.When the crack approaches the edge of the model at a sharp angle, the area of collecting integration points sticks out beyond the model completely on the left side of the crack, as shown in Figure 6a.Theoretically, the crack should continue at the same angle θ as before (about −30 • ), but, as can be seen in Figure 6b, this angle is outside the range of the integration points taken into account.This leads to the fact that fitting any polynomial is highly inaccurate.In such a case, most often, the simulation from this moment gives a very unlikely crack path.
Fortunately, this problem is not so significant, as this phenomenon occurs very close to the edge of the model, i.e., at the very end of the simulation.There is also another problem that seems to have no simple solution.When the crack approaches the edge of the model at a sharp angle, the area of collecting integration points sticks out beyond the model completely on the left side of the crack, as shown in Figure 6a.Theoretically, the crack should continue at the same angle θ as before (about −30°), but, as can be seen in Figure 6b, this angle is outside the range of the integration points taken into account.This leads to the fact that fitting any polynomial is highly inaccurate.In such a case, most often, the simulation from this moment gives a very unlikely crack path.Fortunately, this problem is not so significant, as this phenomenon occurs very close to the edge of the model, i.e., at the very end of the simulation.There is also another problem that seems to have no simple solution.When the crack approaches the edge of the model at a sharp angle, the area of collecting integration points sticks out beyond the model completely on the left side of the crack, as shown in Figure 6a.Theoretically, the crack should continue at the same angle θ as before (about −30°), but, as can be seen in Figure 6b, this angle is outside the range of the integration points taken into account.This leads to the fact that fitting any polynomial is highly inaccurate.In such a case, most often, the simulation from this moment gives a very unlikely crack path.Fortunately, this problem is not so significant, as this phenomenon occurs very close to the edge of the model, i.e., at the very end of the simulation.There is also another problem that seems to have no simple solution.When the crack approaches the edge of the model at a sharp angle, the area of collecting integration points sticks out beyond the model completely on the left side of the crack, as shown in Figure 6a.Theoretically, the crack should continue at the same angle θ as before (about −30°), but, as can be seen in Figure 6b, this angle is outside the range of the integration points taken into account.This leads to the fact that fitting any polynomial is highly inaccurate.In such a case, most often, the simulation from this moment gives a very unlikely crack path.Fortunately, this problem is not so significant, as this phenomenon occurs very close to the edge of the model, i.e., at the very end of the simulation.

Rankine Criterion
Five stress-based material failure criteria were programmed [17][18][19].The first criterion was the Rankine criterion, which was programmed by the authors earlier and has already been used for the authors' previous research.
The criterion envelope in the 2D stress state in the principal stress plane is shown in Figure 7, where σ 1 and σ 3 are the maximum and minimum principal stresses, and f t and f c are the tensile and compressive strength of the material, respectively.For the purposes of the analyses, tensile stresses were assumed to be positive, and compressive stresses were assumed to be negative.In theory, material fracture is influenced by both tensile and compressive strength.However, it was noticed in the simulations that it is not possible for the compressive stresses to exceed the compressive strength before the tensile stresses exceed the tensile strength.Therefore, it was decided to leave the value of material effort µ in the simplest possible form:

Rankine Criterion
Five stress-based material failure criteria were programmed [17][18][19].The first criterion was the Rankine criterion, which was programmed by the authors earlier and has already been used for the authors' previous research.
The criterion envelope in the 2D stress state in the principal stress plane is shown in Figure 7, where σ1 and σ3 are the maximum and minimum principal stresses, and ft and fc are the tensile and compressive strength of the material, respectively.For the purposes of the analyses, tensile stresses were assumed to be positive, and compressive stresses were assumed to be negative.In theory, material fracture is influenced by both tensile and compressive strength.However, it was noticed in the simulations that it is not possible for the compressive stresses to exceed the compressive strength before the tensile stresses exceed the tensile strength.Therefore, it was decided to leave the value of material effort µ in the simplest possible form: It should be remembered that cracking occurs when µ is greater than 1, that is, when the maximum principal stress point lies outside the criterion envelope (σ1 > ft).It is also worth noting that the criterion in this form does not require the compressive strength value.

Coulomb-Mohr Criterion
In the case of the Coulomb-Mohr criterion in the compression-tension quadrant, the envelope in the 2D stress state in the plane of principal stresses takes the shape of a linear function (Figure 8a).The material effort is the ratio of the resultant stress values to the proportional stresses located on the criterion envelope: The intersection of the envelope with the line connecting the center of the coordinate system with the stress data point can be determined using simple formulas: where: η-ratio of compressive strength to tensile strength fc/ft; κ-ratio of maximum principal stresses to minimum σ1/σ3.It should be remembered that cracking occurs when µ is greater than 1, that is, when the maximum principal stress point lies outside the criterion envelope (σ 1 > f t ).It is also worth noting that the criterion in this form does not require the compressive strength value.

Coulomb-Mohr Criterion
In the case of the Coulomb-Mohr criterion in the compression-tension quadrant, the envelope in the 2D stress state in the plane of principal stresses takes the shape of a linear function (Figure 8a).The material effort is the ratio of the resultant stress values to the proportional stresses located on the criterion envelope: The purpose of the discussed method for predicting the direction of crack propagation is also the ability to determine the direction in tasks with a triaxial stress state (for example, axially symmetric tasks), i.e., when the middle principal stresses are different from 0. As can be seen in Figure 8b, a cross-section of the envelope with the σ1, σ3 plane will be different depending on the value of σ2 (an example cross-section for σ2 > 0 is marked with a blue line in this figure), so the method described in the above formulas will be incorrect.However, the material strength can be determined from the following formula: where The intersection of the envelope with the line connecting the center of the coordinate system with the stress data point can be determined using simple formulas: where: η-ratio of compressive strength to tensile strength f c /f t ; κ-ratio of maximum principal stresses to minimum σ 1 /σ 3 .
The purpose of the discussed method for predicting the direction of crack propagation is also the ability to determine the direction in tasks with a triaxial stress state (for example, axially symmetric tasks), i.e., when the middle principal stresses are different from 0. As can be seen in Figure 8b, a cross-section of the envelope with the σ 1 , σ 3 plane will be different depending on the value of σ 2 (an example cross-section for σ 2 > 0 is marked with a blue line in this figure), so the method described in the above formulas will be incorrect.However, the material strength can be determined from the following formula: where The advantage of the above formula is that it works for any stress configuration (compression-tension and tension-tension quadrant) and for a three-dimensional stress state.This formula, due to its relative ease of application to that described by Equation ( 2), was used in the algorithm implemented in the Abaqus ® system.

Drucker-Prager Criterion
The implementation of the Drucker-Prager criterion was almost identical to that described above for the Coulomb-Mohr criterion.Also, in this case, Formula (2) was applied for the 2D stress state, but the coordinates of the intersection of the stress line and the envelope are calculated differently since the shape of the envelope in this coordinate system is an ellipse (Figure 9a):

Ottosen-Podgórski Criterion
The Ottosen-Podgórski criterion [20] was proposed in a form containing three stress tensor invariants: where P(J) is a function describing the cross-section of the criterion in a deviatoric plane where σ0 = const.: where , , C0, C1, C2-parameters obtained from material constants.More information about the discussed criterion can be found in publications [15,20,21] and Appendix A. The shape of the criterion envelope in the three-dimensional stress state is shown in Figure 10a.The material effort itself can be determined from the following ratio: Also, in this case, the envelope differ if the middle principal stresses are different from 0 (Figure 9b).In the case of the Drucker-Prager criterion, a formula analogous to (4) are used: For use in the Abaqus ® system, only the method described by the above formula was programmed for the Drucker-Prager criterion.

Ottosen-Podgórski Criterion
The Ottosen-Podgórski criterion [20] was proposed in a form containing three stress tensor invariants: where P(J) is a function describing the cross-section of the criterion in a deviatoric plane where σ 0 = const.: where ξ, φ, C 0 , C 1 , C 2 -parameters obtained from material constants.More information about the discussed criterion can be found in publications [15,20,21] and Appendix A. The shape of the criterion envelope in the three-dimensional stress state is shown in Figure 10a.
The material effort itself can be determined from the following ratio: where σ 0 , τ 0 -normal and tangential stresses occurring at a given point; σ f , τ f -normal and tangential stresses proportional to the above.This is a point located on the criterion surface.The above stress values are also shown in Figure 10b.

Ottosen-Podgórski Criterion
The Ottosen-Podgórski criterion [20] was proposed in a form containing three stress tensor invariants: where P(J) is a function describing the cross-section of the criterion in a deviatoric plane where σ0 = const.: where , , C0, C1, C2-parameters obtained from material constants.More information about the discussed criterion can be found in publications [15,20,21] and Appendix A. The shape of the criterion envelope in the three-dimensional stress state is shown in Figure 10a.The material effort itself can be determined from the following ratio: where σ0, τ0-normal and tangential stresses occurring at a given point; σf, τf-normal and tangential stresses proportional to the above.This is a point located on the criterion surface.The above stress values are also shown in Figure 10b.

Hoek-Brown Criterion
The Hoek-Brown criterion was originally proposed in 1980 [22] in the form: where A and B are material constants: As for the Coulomb-Mohr and Drucker-Prager criteria, the Hoek-Brown criterion is described in a 2D stress state in the plane of principal stresses [23,24].The shape of the criterion envelope in this plane is shown in Figure 11.Here, too, the material effort could Materials 2024, 17, 3930 9 of 18 be determined using Formula (2), but the point of intersection of the stress line with the envelope could be determined from the formulas obtained after transforming Formula (11): where A and B are material constants: As for the Coulomb-Mohr and Drucker-Prager criteria, the Hoek-Brown criterion is described in a 2D stress state in the plane of principal stresses [23,24].The shape of the criterion envelope in this plane is shown in Figure 11.Here, too, the material effort could be determined using Formula ( 2), but the point of intersection of the stress line with the envelope could be determined from the formulas obtained after transforming Formula (11): Unfortunately, due to insufficient information in the literature, the analysis for the three-dimensional stress state was not performed, so in simulations with three stress directions, the above failure criterion was used for the 2D stress state.Unfortunately, due to insufficient information in the literature, the analysis for the three-dimensional stress state was not performed, so in simulations with three stress directions, the above failure criterion was used for the 2D stress state.

Comparison of All Criteria
Figure 12 summarizes all the analyzed material failure criteria in the plane of principal stresses for the ratio f c /f t = 10.
where A and B are material constants: As for the Coulomb-Mohr and Drucker-Prager criteria, the Hoek-Brown criterion is described in a 2D stress state in the plane of principal stresses [23,24].The shape of the criterion envelope in this plane is shown in Figure 11.Here, too, the material effort could be determined using Formula ( 2), but the point of intersection of the stress line with the envelope could be determined from the formulas obtained after transforming Formula (11): Unfortunately, due to insufficient information in the literature, the analysis for the three-dimensional stress state was not performed, so in simulations with three stress directions, the above failure criterion was used for the 2D stress state.As can be observed, in the compression-tension quadrant, the criteria are slightly different from each other.The Rankine criterion is the upper limit of the envelope, and the Coulomb-Mohr criterion is the lower limit.Hence, the Rankine criterion should give the lowest value of material effort and the Coulomb-Mohr criterion the highest.In the tension-tension quadrant, only the Drucker-Prager criterion differs from the others.For this criterion, a much higher value of the material effort should be obtained than for the other criteria.

Comparison of All Criteria
In order to verify the correctness of the above-described formulas for all criteria, it was decided to determine the material effort values for three real exemplary stress tensors: The following material parameters were selected: f c = 92.56MPa, f t = 3.11 MPa.These are the material parameters of the "Brenna" sandstone, which was the subject of laboratory tests in the authors' previous works.The first stress tensor is a simple example of principal stresses in the tension-compression quadrant.The remaining two are actual stress values read randomly from simulations performed for the purposes of this publication, respectively, for the 2D and the axisymmetric tasks.
Table 1 lists the above-described values of the material effort.In the case of a plane stress state, both ways of determining the coefficient (as for a plane stress state and for a three-dimensional stress state) were the same for the Coulomb-Mohr and Drucker-Prager criteria.However, for tensor C, i.e., for a three-dimensional task, the value of the coefficient determined in the same way as for a flat task was obviously incorrect, but the difference was not very significant.In each case, the material effort obtained the minimum value for the Rankine criterion and the maximum value for the Coulomb-Mohr criterion.For tensors B and C (biaxial stretching), the coefficient values obtained for the Drucker-Prager criterion differed significantly from the other criteria.All the assumptions deduced at the beginning of this chapter are, therefore, true.It may also be interesting to analyze the material effort and, thus, the crack propagation angle based on theoretical data.A Griffith's crack was used, i.e., a crack placed on an infinite disk loaded with a uniform load.A diagram of this theoretical task is shown in Figure 13a.This crack can be loaded perpendicularly-in Mode I of loading (Figure 13b) or longitudinally-in Mode II (Figure 13c).The stresses around the crack tip in both modes were solved by Westergaard [16] and are described by the following formulas:  The stresses around the crack tip in both modes were solved by Westergaard [16] and are described by the following formulas: where K I and K II are the critical stress intensity factors for the material, taken as 1, r-distance from the crack tip, taken as 1, θ-angle around the crack tip.
The principal stresses and octahedral stresses were determined for the above functions, and on their basis, formulas for the material effort functions were determined for all the discussed criteria.The ratio f c /f t was also assumed to be 10.The plots of these functions for both loading modes are shown in Figure 14

Three-Point Bending Test of a Notched Beam
A three-point bending test of a notched beam was the subject of the authors' research.The present work contains a description of laboratory tests of material parameters and the results of the described test.The test scheme is shown in Figure 15a, and the finite element mesh in Figure 15b

Three-Point Bending Test of a Notched Beam
A three-point bending test of a notched beam was the subject of the authors' research.The present work contains a description of laboratory tests of material parameters and the results of the described test.The test scheme is shown in Figure 15a, and the finite element mesh in Figure 15b In computer simulations, the load was controlled by a vertical displacement at a node at the center of the top edge.The sample was modeled in a plane stress state, with supports modeled in the lower left and right corners.The mesh in the crack area was modeled as uniform and rectangular.The size of finite elements ranges from 1 to 15 mm.The obtained force-displacement plots and crack paths are shown in Figure 16.In computer simulations, the load was controlled by a vertical displacement at a node at the center of the top edge.The sample was modeled in a plane stress state, with supports modeled in the lower left and right corners.The mesh in the crack area was modeled as uniform and rectangular.The size of finite elements ranges from 1 to 15 mm.The obtained force-displacement plots and crack paths are shown in Figure 16.
A three-point bending test of a notched beam was the subject of the authors' research.The present work contains a description of laboratory tests of material parameters and the results of the described test.The test scheme is shown in Figure 15a, and the finite element mesh in Figure 15b In computer simulations, the load was controlled by a vertical displacement at a node at the center of the top edge.The sample was modeled in a plane stress state, with supports modeled in the lower left and right corners.The mesh in the crack area was modeled as uniform and rectangular.The size of finite elements ranges from 1 to 15 mm.The obtained force-displacement plots and crack paths are shown in Figure 16.The nature of the damage was the same for all criteria.The maximum force was similar, and only for the default Abaqus ® criterion was overestimated.The Drucker-Prager criterion gave an underestimated result since this was a 2D task, and the stresses around the crack tip were mainly tensile stresses.The crack path for the default Abaqus ® method was only initially consistent with the trajectory predicted by the analytical solutions.After crossing about 2/3 of the path, it unexpectedly turned (Figure 16b), which is probably related to the highly simplified method of estimating the stress field around the crack tip.This change in trajectory led to an increase in the value of the force required for further crack propagation (Figure 16a) and the inability to drive the crack to the top edge of the model.All programmed criteria generated very correct crack paths, maintaining symmetry and brought to the edge of the model.This crack path was consistent with the results obtained in the laboratory, as described in the publication [8].Unfortunately, the maximum force obtained was far from the laboratory test result due to the heterogeneity of the sandstone, as mentioned in a previous publication [25].

Pull-Out Test of a Self-Undercut Anchor
The self-undercutting anchor pull-out test has also been described in previous publications [8].This is an axisymmetric task.The task diagram and the finite element mesh are presented in Figure 17.The material used was the same as for the study above, and the anchor length was h 0 = 10 cm. the heterogeneity of the sandstone, as mentioned in a previous publication [25].

Pull-Out Test of a Self-Undercut Anchor
The self-undercutting anchor pull-out test has also been described in previous publications [8].This is an axisymmetric task.The task diagram and the finite element mesh are presented in Figure 17.The material used was the same as for the study above, and the anchor length was h0 = 10 cm.In the simulations, the load was controlled by vertical displacement.In situ tests were performed on a mass of rock, so the model was chosen to be large enough so the bottom and right edges of the model would not have a significant impact on the fracture.Boundary conditions were created on the lower and right edges of the model, and on the left edge, there was an axis of rotation.The size of the finite element mesh ranged from 2 mm near the expected beginning of the crack to 25 mm.The results of computer simulations, such as the plot of the force-displacement relation and the obtained crack paths, are shown in Figure 18.
The relationship between the pull-out force and the vertical displacement was again very similar for all criteria.Again, the default method gave slightly overestimated results.As for the crack line, again, the default method gave an incorrect result.However, the programmed criteria led the line correctly and very similarly to each other.As proven in previous publications, the crack line and the destructive force for the authors' method coincided with the results of in situ tests [8] (the pull-out force obtained on in situ test was 123 kN).In the simulations, the load was controlled by vertical displacement.In situ tests were performed on a mass of rock, so the model was chosen to be large enough so the bottom and right edges of the model would not have a significant impact on the fracture.Boundary conditions were created on the lower and right edges of the model, and on the left edge, there was an axis of rotation.The size of the finite element mesh ranged from 2 mm near the expected beginning of the crack to 25 mm.The results of computer simulations, such as the plot of the force-displacement relation and the obtained crack paths, are shown in Figure 18.

Pull-Out Test on Rod with Thread
Another test, also described in the authors' previous publication, is a pull-out test of a rod embedded in a concrete cylinder [25].This study was different from the above.Here, there was contact between the rod thread and the concrete, whereas above, the load was carried out only by the pressure of the anchor undercut on the rock.This time, the material parameters were Young's modulus E = 38.000GPa, Poisson's ratio ν = 0.2, compressive strength fc = 44.013MPa, tensile strength ft = 3.300 MPa, critical fracture energy in Mode I GI0 = 0.015 N/mm.The task scheme is shown in Figure 19a.The sample dimensions were width d = 150 mm, height h = 300 mm, and screw mounting depth h0 = 5 cm.A nut was attached to the end of the screw.The load was controlled by vertical displacement through the contact of the thread with the concrete and the pressure of the nut on the concrete.The test was modeled as a 1/4 cross-section of a cylinder, as an axially symmetrical test.On the left and bottom edges, there were axes of symmetry, with the left one being the axis of rotation.The finite element mesh is shown in Figure 19b.The size of finite elements was from 2 to 5 mm.The relationship between the pull-out force and the vertical displacement was again very similar for all criteria.Again, the default method gave slightly overestimated results.As for the crack line, again, the default method gave an incorrect result.However, the programmed criteria led the line correctly and very similarly to each other.As proven in previous publications, the crack line and the destructive force for the authors' method coincided with the results of in situ tests [8] (the pull-out force obtained on in situ test was 123 kN).

Pull-Out Test on Rod with Thread
Another test, also described in the authors' previous publication, is a pull-out test of a rod embedded in a concrete cylinder [25].This study was different from the above.Here, there was contact between the rod thread and the concrete, whereas above, the load was carried out only by the pressure of the anchor undercut on the rock.This time, the material parameters were Young's modulus E = 38.000GPa, Poisson's ratio ν = 0.2, compressive strength f c = 44.013MPa, tensile strength f t = 3.300 MPa, critical fracture energy in Mode I G I0 = 0.015 N/mm.The task scheme is shown in Figure 19a.The sample dimensions were width d = 150 mm, height h = 300 mm, and screw mounting depth h 0 = 5 cm.A nut was attached to the end of the screw.The load was controlled by vertical displacement through the contact of the thread with the concrete and the pressure of the nut on the concrete.The test was modeled as a 1/4 cross-section of a cylinder, as an axially symmetrical test.On the left and bottom edges, there were axes of symmetry, with the left one being the axis of rotation.The finite element mesh is shown in Figure 19b.The size of finite elements was from 2 to 5 mm. a rod embedded in a concrete cylinder [25].This study was different from the above.Here, there was contact between the rod thread and the concrete, whereas above, the load was carried out only by the pressure of the anchor undercut on the rock.This time, the material parameters were Young's modulus E = 38.000GPa, Poisson's ratio ν = 0.2, compressive strength fc = 44.013MPa, tensile strength ft = 3.300 MPa, critical fracture energy in Mode I GI0 = 0.015 N/mm.The task scheme is shown in Figure 19a.The sample dimensions were width d = 150 mm, height h = 300 mm, and screw mounting depth h0 = 5 cm.A nut was attached to the end of the screw.The load was controlled by vertical displacement through the contact of the thread with the concrete and the pressure of the nut on the concrete.The test was modeled as a 1/4 cross-section of a cylinder, as an axially symmetrical test.On the left and bottom edges, there were axes of symmetry, with the left one being the axis of rotation.The finite element mesh is shown in Figure 19b.The size of finite elements was from 2 to 5 mm.The results of computer simulations are shown in Figure 20, which are a forcedisplacement relationship plot and the crack paths.Again, the results of the maximum force and the nature of damage for all criteria were similar, and the fracture lines for all The results of computer simulations are shown in Figure 20, which are a force-displacement relationship plot and the crack paths.Again, the results of the maximum force and the nature of damage for all criteria were similar, and the fracture lines for all programmed criteria were almost identical, while the result for the simulation with the default method differed from the others.Unlike the previous two simulations, the default method brought the crack to the edge of the model.The crack path was similar to the line obtained in the laboratory tests (compare Figures 19a and 20b), and the maximum force was also consistent (Figure 20a).Unfortunately, there was a large difference in the nature of failure (force-displacement relationship).This problem was described in the previous publication [25] and was not caused directly by the programmed criteria.The most likely reason for this was incorrectly assumed fracture energy.For a lower value of the fracture energy in the simulations, the value of the force needed to initiate the fracture decreased.Boundary conditions also turned out to be a problem.Another reason for some inconsistencies was the inability to properly model the connection between the thread and the concrete.The different methods of modeling the load transfer to the concrete did not significantly affect the results, but only in the case of simulations with an elastic rod did the shape of the force-displacement relationship curve turn out to be similar to that obtained in tests.programmed criteria were almost identical, while the result for the simulation with the default method differed from the others.Unlike the previous two simulations, the default method brought the crack to the edge of the model.The crack path was similar to the line obtained in the laboratory tests (compare Figure 19a and Figure 20b), and the maximum force was also consistent (Figure 20a).Unfortunately, there was a large difference in the nature of failure (force-displacement relationship).This problem was described in the previous publication [25] and was not caused directly by the programmed criteria.The most likely reason for this was incorrectly assumed fracture energy.For a lower value of the fracture energy in the simulations, the value of the force needed to initiate the fracture decreased.Boundary conditions also turned out to be a problem.Another reason for some inconsistencies was the inability to properly model the connection between the thread and the concrete.The different methods of modeling the load transfer to the concrete did not significantly affect the results, but only in the case of simulations with an elastic rod did the shape of the force-displacement relationship curve turn out to be similar to that obtained in tests.

Discussion
This paper presents the results of numerical simulations of several laboratory tests performed on specimens made of rock-like materials.These included tests of three-point bending of a notched beam cut from sandstone, the pull-out of a self-undercutting anchor fixed in sandstone, and a pull-out test of a bar embedded in a cylindrical concrete specimen.The simulations were performed using the commercial Simulia Abaqus ®

Conclusions
The numerical simulations presented in this paper, performed with the commercial Abaqus ® system, showed the inadequacies of the simple method of determining the direction of crack propagation modeled by the X-FEM method.This was particularly evident in areas with complex stresses, where I-Mode cannot be used in fracture analysis.In areas with a dominant tensile principal stress, the standard procedure gave satisfactory results.The UDMGINI procedure programmed by the authors, which is responsible for determining the direction of crack propagation, performed much better, allowing a more accurate determination of the crack path, as well as the maximum value of the force destroying the specimen.Simulation of the pull-out test for a sandstone fragment using a self-undercutting anchor clearly showed the advantage of the improved procedure, which takes into account the criterion of the minimum gradient of the material effort field around the crack tip.
Simulations of highly inhomogeneous materials using the isotropic material model available in the standard X-FEM implementation cannot show the exact geometry of the crack path.The apparent effect of scale, which can be characterized by the ratio of the maximum inhomogeneity size to the specimen size (in concrete, it is the maximum aggregate granularity to the minimum specimen size), should be taken into account when creating an FEM model.As it currently stands, the proposed procedure allows for the analysis of cracking in 2D models (including axisymmetric models).It would be desirable to extend it with the possibility of analyzing 3D phenomena.
The problem, as it seems at present, is also the analysis of bifurcation cracks, which is related to the implementation of the X-FEM method available in Abaqus ® , which requires the use of another specialized system or own software.

Figure 1 .
Figure 1.Rotation of stress tensors to principal stresses at integration points in an enriched finite element.

Figure 1 .
Figure 1.Rotation of stress tensors to principal stresses at integration points in an enriched finite element.

Figure 2 .
Figure 2. Dependence between crack propagation angle α and initial crack angle β for a Griffith's crack for programmed criteria and the maximum energy release rate (MERR) criterion.O-P: minimum gradient of the effort function for the Ottosen-Podgórski failure criterion, maxps: maximum principal stress, α = β: direction normal to applied stress.

Figure 2 .
Figure 2. Dependence between crack propagation angle α and initial crack angle β for a Griffith's crack for programmed criteria and the maximum energy release rate (MERR) criterion.O-P: minimum gradient of the effort function for the Ottosen-Podgórski failure criterion, maxps: maximum principal stress, α = β: direction normal to applied stress.

Figure 3 .
Figure 3. Method of determining the direction of crack propagation: (a) example of selecting integration points taken into account in determining the material effort plot; (b) example of the relationship between material effort and the angle around the crack tip.

Figure 3 .
Figure 3. Method of determining the direction of crack propagation: (a) example of selecting integration points taken into account in determining the material effort plot; (b) example of the relationship between material effort and the angle around the crack tip.

Figure 4 .Figure 5 .
Figure 4. Example of the value of material effort around the crack tip in the case of a crack perpendicular to the edge of the model: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 6 .
Figure 6.Example of the value of material effort around the crack tip in the case of a crack at a high angle to the edge of the model and with a crack tip close to the edge: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 4 .
Figure 4. Example of the value of material effort around the crack tip in the case of a crack perpendicular to the edge of the model: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 4 .Figure 5 .
Figure 4. Example of the value of material effort around the crack tip in the case of a crack perpendicular to the edge of the model: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 6 .
Figure 6.Example of the value of material effort around the crack tip in the case of a crack at a high angle to the edge of the model and with a crack tip close to the edge: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 5 .
Figure 5. Example of the value of material effort around the crack tip in the case of a crack perpendicular to the edge of the model and with the crack tip close to the edge: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 4 .Figure 5 .
Figure 4. Example of the value of material effort around the crack tip in the case of a crack perpendicular to the edge of the model: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 6 .
Figure 6.Example of the value of material effort around the crack tip in the case of a crack at a high angle to the edge of the model and with a crack tip close to the edge: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 6 .
Figure 6.Example of the value of material effort around the crack tip in the case of a crack at a high angle to the edge of the model and with a crack tip close to the edge: (a) area of collecting the integration points; (b) relationship between material effort and the angle around the crack tip.

Figure 7 .
Figure 7. Envelope of the Rankine criterion in the plane of principal stresses.

Figure 7 .
Figure 7. Envelope of the Rankine criterion in the plane of principal stresses.

Figure 8 .
Figure 8. Envelope of the Coulomb-Mohr criterion in the plane of principal stresses: (a) in a twodimensional stress state; (b) in a three-dimensional stress state.

Figure 8 .
Figure 8. Envelope of the Coulomb-Mohr criterion in the plane of principal stresses: (a) in a twodimensional stress state; (b) in a three-dimensional stress state.

Figure 9 .
Figure 9. Envelope of the Drucker-Prager criterion in the plane of principal stresses: (a) in a twodimensional stress state; (b) in a three-dimensional stress state.

Figure 9 .
Figure 9. Envelope of the Drucker-Prager criterion in the plane of principal stresses: (a) in a twodimensional stress state; (b) in a three-dimensional stress state.

Figure 9 .
Figure 9. Envelope of the Drucker-Prager criterion in the plane of principal stresses: (a) in a twodimensional stress state; (b) in a three-dimensional stress state.

Figure 11 .
Figure 11.The envelope of the Hoek-Brown criterion in the plane of principal stresses in a plane stress state.

Figure 12
Figure 12 summarizes all the analyzed material failure criteria in the plane of principal stresses for the ratio fc/ft = 10.

Figure 12 .
Figure 12.Envelopes of all failure criteria in the plane of principal stresses.

Figure 11 .
Figure 11.The envelope of the Hoek-Brown criterion in the plane of principal stresses in a plane stress state.

Figure 11 .
Figure 11.The envelope of the Hoek-Brown criterion in the plane of principal stresses in a plane stress state.

Figure 12
Figure 12 summarizes all the analyzed material failure criteria in the plane of principal stresses for the ratio fc/ft = 10.

Figure 12 .
Figure 12.Envelopes of all failure criteria in the plane of principal stresses.

Figure 12 .
Figure 12.Envelopes of all failure criteria in the plane of principal stresses.

Figure 13 .
Figure 13.The crack analyzed by Griffith (a), and crack loading schemes Mode I (b), Mode II (c).

Figure 13 .
Figure 13.The crack analyzed by Griffith (a), and crack loading schemes Mode I (b), Mode II (c).

Figure 14 .
Figure 14.Plots of the material effort for stresses around the tip of the Griffith's crack: (a) for Mode I; (b) for Mode II.

Figure 14 .
Figure 14.Plots of the material effort for stresses around the tip of the Griffith's crack: (a) for Mode I; (b) for Mode II.

Figure 15 .
Figure 15.Three-point bending test of a beam with a notch: (a) scheme of the task; (b) finite element mesh.

Figure 15 .
Figure 15.Three-point bending test of a beam with a notch: (a) scheme of the task; (b) finite element mesh.

Figure 15 .
Figure 15.Three-point bending test of a beam with a notch: (a) scheme of the task; (b) finite element mesh.

Figure 16 .
Figure 16.Simulation results of three-point bending of a beam with a notch: (a) force-displacement plot; (b) crack paths.

Figure 17 .
Figure 17.Pull-out test of a self-undercutting anchor: (a) scheme of the task; (b) finite element mesh.

Figure 17 .
Figure 17.Pull-out test of a self-undercutting anchor: (a) scheme of the task; (b) finite element mesh.

Figure 18 .
Figure 18.Simulation results of the pull-out test of a self-undercutting anchor: (a) forcedisplacement plot; (b) crack paths after simulation compared with a rock sample pulled out with a self-undercutting anchor.

Figure 18 .
Figure 18.Simulation results of the pull-out test of a self-undercutting anchor: (a) force-displacement plot; (b) crack paths after simulation compared with a rock sample pulled out with a self-undercutting anchor.

Figure 19 .
Figure 19.Pull-out test of a rod embedded in concrete: (a) scheme of the task; (b) finite element mesh; (c) laboratory test; (d) exemplary sample after laboratory test.

Figure 19 .
Figure 19.Pull-out test of a rod embedded in concrete: (a) scheme of the task; (b) finite element mesh; (c) laboratory test; (d) exemplary sample after laboratory test.

Figure 20 .
Figure 20.Simulation results of the pull-out test of a rod embedded in concrete: (a) forcedisplacement plot; (b) crack paths.

Figure 20 .
Figure 20.Simulation results of the pull-out test of a rod embedded in concrete: (a) force-displacement plot; (b) crack paths.

Table 1 .
Material effort for three exemplary stress tensors for all criteria.